% DETERMINE WHICH TRIALS TO OMIT FROM GRAND AVERAGES

omit = [];
p = 1;
for t = 1:numel(files)
    if trialsummary(1,t).RMSE > outlier  % adjust RMSE cut off if nec
        omit(p) = t;
        p = p + 1;
    end
end

% CALCULATE AVERAGE TRACES & RMSE'S OVER ALL VALID TRIALS (3 CYCLES)

fn = 0;
fs = 0;
sn = 0;
ss = 0;
Empty = zeros(18000,1);

FN.X = Empty;
FS.X = Empty;
SN.X = Empty;
SS.X = Empty;

FN.allerrsq = [];
FS.allerrsq = [];
SN.allerrsq = [];
SS.allerrsq = [];

FN.allerrsqocc = [];
FS.allerrsqocc = [];
SN.allerrsqocc = [];
SS.allerrsqocc = [];

FN.allerrsqfix = [];
FS.allerrsqfix = [];
SN.allerrsqfix = [];
SS.allerrsqfix = [];

for t = 1:numel(files)
    if ismember(t,omit) == 0  % is trial valid?
        switch trialsummary(1,t).speed.deg
            case 22
                if trialsummary(1,t).stimulus.smooth == 0
                    FN.X(:,fn+1) = trialsummary(1,t).data.X;
                    FN.velnosac(:,fn+1) = trialsummary(1,t).data.vel;
                    FN.allerrsq = [FN.allerrsq;...
                        trialsummary(1,t).errsq];
                    FN.allerrsqocc = [FN.allerrsqocc;...
                        trialsummary(1,t).errsqocc];
                    FN.allerrsqfix = [FN.allerrsqfix;...
                        trialsummary(1,t).errsqfix];
                    fn = fn + 1;
                else
                    FS.X(:,fs+1) = trialsummary(1,t).data.X;
                    FS.velnosac(:,fs+1) = trialsummary(1,t).data.vel;
                    FS.allerrsq = [FS.allerrsq;...
                        trialsummary(1,t).errsq];
                    FS.allerrsqocc = [FS.allerrsqocc;...
                        trialsummary(1,t).errsqocc];
                    FS.allerrsqfix = [FS.allerrsqfix;...
                        trialsummary(1,t).errsqfix];
                    fs = fs + 1;
                end
            case 18
                if trialsummary(1,t).stimulus.smooth == 0
                    SN.X(:,sn+1) = trialsummary(1,t).data.X;
                    SN.velnosac(:,sn+1) = trialsummary(1,t).data.vel;
                    SN.allerrsq = [SN.allerrsq;...
                        trialsummary(1,t).errsq];
                    SN.allerrsqocc = [SN.allerrsqocc;...
                        trialsummary(1,t).errsqocc];
                    SN.allerrsqfix = [SN.allerrsqfix;...
                        trialsummary(1,t).errsqfix];
                    sn = sn + 1;
                else
                    SS.X(:,ss+1) = trialsummary(1,t).data.X;
                    SS.velnosac(:,ss+1) = trialsummary(1,t).data.vel;
                    SS.allerrsq = [SS.allerrsq;...
                        trialsummary(1,t).errsq];
                    SS.allerrsqocc = [SS.allerrsqocc;...
                        trialsummary(1,t).errsqocc];
                    SS.allerrsqfix = [SS.allerrsqfix;...
                        trialsummary(1,t).errsqfix];
                    ss = ss + 1;
                end
        end
        
    end
end

FN.mean = nanmean(FN.X,2);
FN.mean = FN.mean(1:12519);
FS.mean = nanmean(FS.X,2);
FS.mean = FS.mean(1:12519);
SN.mean = nanmean(SN.X,2);
SN.mean = SN.mean(1:15300);
SS.mean = nanmean(SS.X,2);
SS.mean = SS.mean(1:15300);

FN.total = fn;
FS.total = fs;
SN.total = sn;
SS.total = ss;

FN.time = 1:12519;
FS.time = 1:12519;
SN.time = 1:15300;
SS.time = 1:15300;

FN.RMSE = sqrt(nanmean(FN.allerrsq));
FS.RMSE = sqrt(nanmean(FS.allerrsq));
SN.RMSE = sqrt(nanmean(SN.allerrsq));
SS.RMSE = sqrt(nanmean(SS.allerrsq));
Stats.RMSE(:,1) = [FN.RMSE, FS.RMSE, SN.RMSE, SS.RMSE];

FN.RMSEocc = sqrt(nanmean(FN.allerrsqocc));
FS.RMSEocc = sqrt(nanmean(FS.allerrsqocc));
SN.RMSEocc = sqrt(nanmean(SN.allerrsqocc));
SS.RMSEocc = sqrt(nanmean(SS.allerrsqocc));
Stats.RMSE(:,2) = [FN.RMSEocc, FS.RMSEocc, SN.RMSEocc, SS.RMSEocc];

FN.RMSEfix = sqrt(nanmean(FN.allerrsqfix));
FS.RMSEfix = sqrt(nanmean(FS.allerrsqfix));
SN.RMSEfix = sqrt(nanmean(SN.allerrsqfix));
SS.RMSEfix = sqrt(nanmean(SS.allerrsqfix));
Stats.RMSE(:,3) = [FN.RMSEfix, FS.RMSEfix, SN.RMSEfix, SS.RMSEfix];

% CALCULATE AVERAGE OVER ALL CYCLES

FN.one(:,1) = FN.mean(1:4173,1); 
FN.one(:,2) = FN.mean(4174:8346,1);
FN.one(:,3) = FN.mean(8347:12519,1);
FN.one = nanmean(FN.one,2);

FS.one(:,1) = FS.mean(1:4173,1); 
FS.one(:,2) = FS.mean(4174:8346,1);
FS.one(:,3) = FS.mean(8347:12519,1);
FS.one = nanmean(FS.one,2);

SN.one(:,1) = SN.mean(1:5100,1); 
SN.one(:,2) = SN.mean(5101:10200,1);
SN.one(:,3) = SN.mean(10201:15300,1);
SN.one = nanmean(SN.one,2);

SS.one(:,1) = SS.mean(1:5100,1); 
SS.one(:,2) = SS.mean(5101:10200,1);
SS.one(:,3) = SS.mean(10201:15300,1);
SS.one = nanmean(SS.one,2);
      
% CALCULATE MEAN VELOCITIES EXCL SACCADES & SMOOTH

FN.meanvelnosac = nanmean(FN.velnosac,2);
FN.meanvelnosac = smooth(FN.meanvelnosac(1:12519),smve);
FS.meanvelnosac = nanmean(FS.velnosac,2);
FS.meanvelnosac = smooth(FS.meanvelnosac(1:12519),smve);
SN.meanvelnosac = nanmean(SN.velnosac,2);
SN.meanvelnosac = smooth(SN.meanvelnosac(1:15300),smve);
SS.meanvelnosac = nanmean(SS.velnosac,2);
SS.meanvelnosac = smooth(SS.meanvelnosac(1:15300),smve);

FN.onevelnosac = (FN.meanvelnosac(1:4173,1) + ...
                  FN.meanvelnosac(4174:8346,1) + ...
                  FN.meanvelnosac(8347:12519,1))./3;
SN.onevelnosac = (SN.meanvelnosac(1:5100,1) + ...
                  SN.meanvelnosac(5101:10200,1) + ...
                  SN.meanvelnosac(10201:15300,1))./3;
FS.onevelnosac = (FS.meanvelnosac(1:4173,1) + ...
                  FS.meanvelnosac(4174:8346,1) + ...
                  FS.meanvelnosac(8347:12519,1))./3;
SS.onevelnosac = (SS.meanvelnosac(1:5100,1) + ...
                  SS.meanvelnosac(5101:10200,1) + ...
                  SS.meanvelnosac(10201:15300,1))./3;

% CALCULATE MEAN VELOCITIES INCL SACCADES (smoothing sm = 100)

FN.meanvel = velocity(FN.mean,smvi,scnw,scnwcm,view_dist);
SN.meanvel = velocity(SN.mean,smvi,scnw,scnwcm,view_dist);
FS.meanvel = velocity(FS.mean,smvi,scnw,scnwcm,view_dist);
SS.meanvel = velocity(SS.mean,smvi,scnw,scnwcm,view_dist);

FN.meanvel = FN.meanvel(1:12519);
FS.meanvel = FS.meanvel(1:12519);
SN.meanvel = SN.meanvel(1:15300);
SS.meanvel = SS.meanvel(1:15300);

FN.onevel = velocity(FN.one,smvi,scnw,scnwcm,view_dist);
SN.onevel = velocity(SN.one,smvi,scnw,scnwcm,view_dist);
FS.onevel = velocity(FS.one,smvi,scnw,scnwcm,view_dist);
SS.onevel = velocity(SS.one,smvi,scnw,scnwcm,view_dist);

FN.onevel = FN.onevel(1:4173,1);
FS.onevel = FS.onevel(1:4173,1);
SN.onevel = SN.onevel(1:5100,1);
SS.onevel = SS.onevel(1:5100,1);

% CALCULATE LAGS

FN.lag(1:2086,1) = target22.X(1:2086,1) - FN.one(1:2086,1);
FN.lag(2087:4173,1) = FN.one(2087:4173,1) - target22.X(2087:4173,1);
FS.lag(1:2086,1) = target22.X(1:2086,1) - FS.one(1:2086,1);
FS.lag(2087:4173,1) = FS.one(2087:4173,1) - target22.X(2087:4173,1);

SN.lag(1:2550,1) = target18.X(1:2550,1) - SN.one(1:2550,1);
SN.lag(2551:5100,1) = SN.one(2551:5100,1) - target18.X(2551:5100,1);
SS.lag(1:2550,1) = target18.X(1:2550,1) - SS.one(1:2550,1);
SS.lag(2551:5100,1) = SS.one(2551:5100,1) - target18.X(2551:5100,1);